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Abstract. A quantum shuttle is an archetypical nanoelectromechanical device, 
where the mechanical degree of freedom is quantized. Using a full-scale numerical 
solution of the generalized master equation describing the shuttle, we have recently 
shown [Novotny et al, Phys. Rev. Lett. 92, 248302 (2004)] that for certain limits 
of the shuttle parameters one can distinguish three distinct charge transport 
mechanisms: (i) an incoherent tunneling regime, (ii) a shuttling regime, where 
the charge transport is synchronous with the mechanical motion, and (iii) a 
coexistence regime, where the device switches between the tunneling and shuttling 
regimes. While a study of the cross-over between these three regimes requires the 
full numerics, we show here that by identifying the appropriate time-scales it 
is possible to derive vastly simpler equations for each of the three regimes. The 
simplified equations allow a clear physical interpretation, are easily solved, and are 
in good agreement with the full numerics in their respective domains of validity. 



1. Introduction 

A generic example of a nanoelectromechanical (NEMS) device is given by the charge 
shuttle (originally proposed by Gorelik et al. P): a movable single-electron device, 
working in the Coulomb blockade regime, which can exhibit regular charge transport, 
where one electron within each mechanical oscillation cycle is transported from the 
source to the drain - see Figure 1 for a schematic illustration (see also Ref. 2 , which 
contains an illustrative computer animation). 

The device shown in Fig. ^ can exhibit a number of different charge transport 
mechanisms, to be discussed below, and transitions between the various regimes can be 
induced by varying a control parameter, such as the applied bias, or the mechanical 
damping. Since its introduction, the charge shuttle has inspired a large number of 
theoretical papers (see, e.g., Refs. 011001301111113 IIIl)- To the best of our 
knowledge, a clear-cut experimental demonstration of a shuttling transition is not yet 
available, though significant progress has been made, such as the driven shuttle of 
Erbe et al. or the nanopillars studied by Scheible and Blick In the present 
paper, we study the quantum shuttle, i.e. a device where also the mechanical motion 
needs to be quantized (the physical condition for this to happen is A ~ Xzp, where 
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Figure 1. Schematic representation of the single-dot shuttle: electrons 
tunnel from the left lead at chemical potential (hl) to the quantum dot and 
eventually to the right lead at lower chemical potential hr. The position 
dependent tunneling amplitudes are indicated. X is the displacement from 
the equilibrium position. The springs represent the harmonic potential in 
which the central dot can move. 



A is the tunneling length describing the exponential decay of the wave-functions into 
vacuum, and Xzp is the quantum-mechanical zero-point amplitude). In contrast to 
many of the earlier theoretical papers on quantum shuttles [HI EI El, where extensive 
numerical calculations were employed, the main aim here is to develop simplified 
models which allow significant analytic progress, and hence lead to a transparent 
physical interpretation. From the previous numerical studies we know that there 
are (at least) three well-defined transport regimes for a quantum nanomechanical 
device: (i) the tunneling regime, (ii) the shuttle regime, and (iii) the coexistence 
regime. As we show in subsequent sections, each of these regimes is characterized by 
certain inequalities governing the various time-scales, and a systematic exploitation 
of these inequalities allows us then to develop the aforementioned simplified models. 
In all three cases we will compare the predictions of the simplified models to the 
ones obtained with the full numerics. While in most cases the comparison is very 
satisfactory, we do not always observe quantitative agreement; these discrepancies are 
analyzed and directions for future work are indicated. 

The paper is organized as follows. In Sections 2, 3 and 4 we briefly introduce the 
microscopic model for the quantum shuttle, introduce the Klein-Kramers equations 
for the Wigner functions, and summarize the phenomenology extracted from previous 
numerical studies, respectively. Section 5 contains the main results of this paper, i.e., 
the derivations and analysis of the simplified models for the three transport regimes. 
We end the paper with a short conclusion of the main results. 

2. The Single-Dot Quantum Shuttle 

The Single-Dot Quantum Shuttle (SDQS) consists of a movable quantum dot (QD) 
suspended between source and drain leads (see Fig.^). The center of mass of the QD 
is confined to a potential that, at least for small displacements from its equilibrium 
position, can be considered harmonic. Due to its small geometric size, the QD has 
a very small capacitance and thus a charging energy that may exceed the thermal 
energy ksT (approaching room temperature in the most recent realizations For 
this reason we assume that only one excess electron can occupy the device (Coulomb 
blockade) and we describe the electronic state of the central dot as a two-level 
system (empty/charged). Electrons can tunnel between leads and dot with tunneling 
amplitudes which depend exponentially on the position of the central island due to the 
decreasing/increasing overlapping of the electronic wave functions. The Hamiltonian 
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of the model reads: 

H = Hsys + -ff leads + -ffbath + -fftun + -f^int (1) 



where 



-ffsys = ^ h imcj^x^ + (ei - e£x)c\c-^ 

2m 2 

k (2) 
Htun = Y^[Ti{x)clci + T^(x)4^ci] + /i.e. 

k 

-ffbath + Hini = gcncric heat bath 

The hat over the position and momentum {x,p) of the dot indicates that they are 
operators since the mechanical degree of freedom is quantized. Using the language of 
quantum optics wc call the movable grain alone the system. This is then coupled to 
two electric baths (the leads) and a generic heat bath. The system is described by a 
single electronic level of energy ei and a harmonic oscillator of mass m and frequency 
Lo. When the dot is charged the electrostatic force (eS) acts on the grain and gives 
the electrical influence on the mechanical dynamics. The electric field £ is generated 
by the voltage drop between left and right lead. In our model, though, it is kept as an 
external parameter, also in view of the fact that wc will always assume the potential 
drop to be much larger than any other energy scale of the system (with the only 
exception of the charging energy of the dot). 

The leads are Fermi seas kept at two different chemical potentials (/i^ and /i^) 
by the external applied voltage (AV^ = (^^ — jiiij/e ) and all the energy levels of the 
system lie well inside the bias window. The oscillator is immersed into a dissipative 
environment that we model as a collection of bosons coupled to the oscillator by a 
weak bilinear interaction: 



-ffbath = ^ hVqdq^dq 



q 



-H'int = %\/ — i— * (C^q + C^q^) 



(3) 



where the operator creates a bath boson with wave number q. The damping rate 
is given by: 

7(u;) = 2Trg^D{i^) (4) 

where D{u}) is the density of states for the bosonic bath at the frequency of the system 
oscillator. A bath that generates a frequency independent 7 is called Ohmic. 

The coupling to the electric baths is introduced with the tunneling Hamiltonian 
iJtun- The tunneling amplitudes Ti{x) and Tr{x) depend exponentially on the position 
operator x and represent the mechanical feedback on the electrical dynamics: 

Ti^r{x) = ti^r exp(=F^/A) (5) 

where A is the tunneling length. The tunneling rates from and to the leads (F^^ij) can 
be expressed in terms of the amplitudes: 

Tl,r = {Tl,r{^)) = /^DL,Rexp (^Ty) \f4,rA (6) 
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where are the densities of states of the left and right lead respectively and the 

average is taken with respect to the quantum state of the oscillator. 

The model has three relevant time scales: the period of the oscillator 27r/(jj, the 
inverse of the damping rate I/7 and the average injection/ejection time I/Tl^h. It 
is possible also to identify three important length scales: the zero point uncertainty 

Xzp = \J the tunneling length A and the displaced oscillator equilibrium position 

d — The ratios between the time scales and the ratios between length scales 

distinguish the different operating regimes of the SDQS. 



3. Klein-Kramers Equation 

The shuttle dynamics has an appealing simple classical interpretation: the name 
"shuttle" suggests the idea of sequential and periodical loading, mechanical transport 
and unloading of electrons between a source and a drain lead. Motivated by the 
possibility of observing signatures of quantum dynamics of the mechanical degree of 
freedom for a nanoscale shuttle, we decided, following the suggestion of Armour and 
MacKinnon 3 , to explore a system with a quantized oscillator. We express our results 
in terms of the Wigner function because in this way we can simultaneously keep the 
intuitive phase-space picture and handle the quantum-classical correspondence 

The phase space of the shuttle device is spanned by the triplet charge - position 
- momentum. Correspondingly, the Wigner function is constructed from the reduced 
density matrix an (i = 0, 1 indicates the empty and charged states respectively): 



Wniq.p.t) 



1 

27r^ 



+ 00 



diiq 



exp 



h 



(7) 



where the reduced density matrix a is defined as the trace over the mechanical and 
thermal baths of the full density matrix: 



(7 = TtbIp} 



(8) 



The dynamics of the shuttle device is then completely described by the equation 
of motion for the Wigner distribution [HI I14|: 



dt 



dWu 
dt 



muj q 



d_ 

dp 



p d d 
m oq op 



•ymhuj ub 



1\ &^ 

^ 2) dp' 

2« Q2n^ 



00 



00 



mijJ'{q — d) 



d 



dp 

TLe-^'i/'Woo - 



p d 
m dq 



d 
dp 



jmhu! Ub 



Qp2n 

1 
2 



dp'' 



(9) 



^ (2n)! 

n=0 ^ ' 



2" Ql-^Wr^ 



dp 



2n 



where {q, p) are the position and momentum coordinates of the mechanical phase space 
and is the Bose distribution calculated at the natural frequency of the harmonic 
oscillator. Only the diagonal charge states enter the Klein-Kramers equations Q: 
the off-diagonal charge states vanish given the incoherence of the leads and arc thus 
excluded from the dynamics. 
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We distinguish in Eqs. © contributions of different physical origin: the coherent 
terms that govern the dynamics of the (shifted) harmonic osciUator, the dissipative 
terms proportional to the mechanical damping constant 7 and, finally the driving 
terms proportional to the bare tunneling rates T^^n. 

The ability of the formalism to treat the quantum-classical correspondence is 
explicit in Eqs. ©: given a length, a mass, and a time scale for the system we can 
rescale the phase space coordinates and an expansion in H/Ssys will appear where Ssys 
is the typical action of the system. Classical systems have a large action S'sys 3> h and 
only the first term (h/Ssys ~^ 0) in the expansion is relevant. In the opposite limit 
S'sys ~ fi the full expansion should be considered. 

4. Phenomenology 

The stationary solution of the Klein-Kramers equations for the Wigner distributions 
(O describes the average long time behavior of the shuttling device. Information about 
the different long time operating regimes can be extracted from the distribution itself 
or from the experimentally accessible stationary current and zero frequency current 
noise. 

The mechanical damping rate 7 is the control parameter of our analysis. At high 
damping rates the total Winger distribution is concentrated around the origin of the 
phase space and represents the harmonic oscillator in its ground state. While reducing 
the mechanical damping a ring develops and, after a short coexistence, the central 
"dot" eventually disappears (Figure ^ . The ring is the noisy representation of the 
low damping limit cycle trajectory (shuttling) that develops from the high damping 
equilibrium position (tunneling). Equilibrium and limit cycle dynamics coexist in 
the intermediate damping bistable configuration where the system randomly switches 
between tunneling and shuttling regimes. The charge resolved Winger distributions 
Woo and Wn also reveal the charge-position (momentum) correlation typical of the 
shuttling regime: for negative displacements and positive momentum (i.e. leaving the 
source lead) the dot is prevalently charged while it is empty for positive displacements 
and negative momentum (coming from the drain lead). 

Also the stationary current and the current noise (expressed in terms of the 
Fano factor) show distinctive features for the different operating regimes. At high 
damping rates the shuttling device behaves essentially like the familiar double-barrier 
system since the dot is (almost) static and far from both electrodes. The current is 
determined essentially by the bare tunneling rates Pl.t? and the Fano factor differs only 
slightly from the values found for resonant tunneling devices {F = 1/2 for a symmetric 
device). At low damping rates the current saturates at one electron per mechanical 
cycle (corresponding to current I/euj = l/2n) since the electrons are shuttled one by 
one from the source to the drain lead by the oscillating dot while the extremely sub- 
poissonian Fano factor reveals the deterministic character of this electron transport 
regime. The fingerprint of the coexistence regime, at intermediate damping rates, 
is a substantial enhancement of the Fano factor. The current interpolates smoothly 
between the shuttling and tunneling limiting values (Figure OJ- 

The cross-over damping rate is determined by the effective tunneling rates of the 
electrons. We get the following physical picture: every time an electron jumps on 
the movable grain the grain is subject to the electrostatic force e£ that accelerates 
it towards the drain lead. Energy is pumped into the mechanical system and the 
dot starts to oscillate. If the damping is high compared to the tunneling rates the 



Simple models for the quantum shuttle 6 



y = 0.025 y = 0.029 Y= 0.1 




Figure 2. Charge resolved Wigner function distributions for different 
mechanical damping rates (horizontal axis: coordinate in units of xq = 
^h/mjjj; vertical axis: momentum in h/xo). The rows represent from top 
to bottom the empty (Woo), charged (Wn) and total (Wtot = Woo + Wn) 
Wigner distribution respectively. The columns represent from left to right 
the shuttling (7 = 0.025 lj), coexistence (7 = 0.029 a;) and tunneling regime 
(7 = 0.1 tj) respectively. The Figure is partially reproduced from 



oscillator dissipates this energy into the environment before the next tunneling event: 
on average the dot remains in its ground state. On the other hand, for very small 
damping the relaxation time of the oscillator is long and multiple "forcing events" 
occur before the relaxation takes place. This continuously drives the oscillator away 
from equilibrium and a stationary state is reached only when the energy pumped per 
cycle into the system is dissipated during the same cycle in the environment. 

5. Simplified models 

We qualitatively described in the previous section three possible operating regimes 
for shuttle-devices. The specific separation of time scales allows us to identify the 
relevant variables and describe each regime by a specific simplified model. Models for 
the tunneling, shuttling and coexistence regime are analyzed separately in the three 
following subsections. We also give a comparison with the full description in terms of 
Wigner distributions, current and currcnt-noisc to illustrate how the models capture 
the relevant dynamics. 
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Figure 3. Left panel - Stationary current for the SDQS vs. damping 7. 
The mechanical dissipation rate 7 and the electrical rate T — Fl = r_R are 
given in units of the mechanical frequency oj, the tunneling length I in terms 
of xq = ^hlmuj. The other parameters are d = 0.5x0 and T — 0. The 
current saturates in the shuttling (low damping) regime to one electron per 
cycle independently from the parameters while is substantially proportional 
to the bare electrical rate T = Tl ~ Fr in the tunneling regime (high 
damping). Right panel - Fano factor for the SDQS vs. damping 7. The 
curves correspond to the same parameters of the left panel. The very low 
noise in the shuttling (low damping) regime is a sign of ordered transport. 
The huge super-poissonian Fano factors correspond to the onset of the 
coexistence regime. The Figure is taken from '7L 



5.1. Renormalized resonant tunneling 

The electrical dynamics has the longest time-scale in the tunneling regime since 
the mechanical relaxation time (which is much longer than the oscillation period) 
is much shorter than the average injection or ejection time. Because of this time-scale 
separation, the observation of the device dynamics would most of the time show two 
mechanically frozen states: 

0. Empty dot in the ground state 

1. Charged dot moved to the shifted equilibrium position by the constant electrostatic 

force e£. 

We combine this observation with a quantum description of the mechanical oscillator 
and possible thermal noise under the assumption that the reduced density matrix of 
the device can be written in the form: 

o-oo(i) =Poo(i)CTth(0) ^^^^ 
(Tii(i) = pii(i)(Tth(e£) 

where 

g-/3(Ho.c-.?'2;) 

is the thermal density matrix of a harmonic oscillator subject to an external force T. 
The functions poo (0 ^'i'^ Pn (0 represent the probability to find the system respectively 
in the state or 1 , respectively. The equations of motion for the probabilities pu (t) 
can be derived by inserting the assumption (|10|l in the definition ^ and taking the 
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integral over the mechanical degrees of freedom in the corresponding Klein-Kramers 
equations lO- This results in the rate equations 

^ f PoQ \ ^ f fflpii -fipoo \ = r[ PoQ \ /12^ 
dt \pii J \ TLPoo-TRPn y - ^ ' ^ ^ 

where 



dt \ Pu J V TlPoo - TrPh J V ^'ii 

Tl = riTrmoch |o-th(0)e"T^| ^Tl J dqdpe^~rWth{q,p) 
fi? = Tr,„cch |crth(ef )eT^| = T/j J dqdpe^Wth{q - rf,p) 



(13) 



are the renormalized injection and ejection rates and Trmech indicates the trace over the 
mechanical degrees of freedom of the device. We have also introduced the Liouvillean 
operator: 

The thermal equilibrium Wigner function Wthiq — d,p) is defined as the Wigner 
representation of the thermal equilibrium density matrix ath{e£) — o'th(™'-^^d): 



q — a 



£muj 



(15) 



where i = \J ^m^^'^^B + 1) reduces to the zero point uncertainty length x^p — \J 
in the zero temperature limit. In the high temperature limit fc^T 3> fiw, tends to 
the thermal length Zth ~ yJksT / {muo'^). Using p5(l in 113|l gives the renormalized 
rates: 

TL^T.eA^y 
tR = TRe^^^iir 

Equations H12|l describe the dynamics of a resonant tunneling device. All the effects 
of the movable grain are contained in the effective rates Tl,Tr. As expected the 
ejection rate is modified by the "classical" shift d of the equilibrium position due to 
the electrostatic force on the charged dot. Note that both rates are also enhanced by 
the fuzziness in the position of the oscillator due to thermal and quantum noise. The 
relevance of this correction is given by the ratio between £ and the tunneling length I. 

5.1.1. Phase-space distribution The phase space distribution for the stationary state 
of the simplified model for the tunneling regime is built on the Wigner representation 
of the thermal density matrix tJth and the stationary solution of the system (|12|l for 
the occupation pa of the electromechanical states i: 

Wg'-'iq.p) = ^^^^WtUq,p) 

^L + Tr ^^^^ 

W!l-'{q,p) = - Wthiq - d,p) 

1 L + 1 fl 

The stationary distribution of the tunneling model is determined by the length £ and 
associated momentum moj£, the equilibrium position shift d, the tunneling length I 
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Figure 4. Comparison between the numerical and the analytical results 
for the Wigner distribution functions. The coordinate (left) or momentum 
(right) cuts always cross the maximum of the distribution. The green (blue) 
circles are numerical results in the charged (empty) dot configuration, and 
the full lines represent the analytical calculations. The parameters are: 
Tl=Vr = 0.01 a;, 7 = 0.25 tj, d = 0.5 xq, / = 2 a;o, T = 0. We also plotted 
with dots the numerical results for T — 0.001 cjJ. 



and the ratio between left and right bare electrical rates Tl/Th. The mechanical 
relaxation rate 7 drops out from the solution and only sets the range of applicability 
of the simplified model. 

In Figures 21 and 13 we compare the Wigner ftmctions calculated both analytically 
and numerically in the tunneling regime. They show in general a good agreement 
(Figure 21). The matching is further improved when reducing the bare injection rate 
Tl.t? = O.OOlw thus enlarging the time scale separation r^ ji <C 7 typical of the 
tunneling regime. The temperature dependence of the stationary Wigner function 
distribution (Figure verifies the scaling given by the temperature dependent length 
L 




Figure 5. Tunneling Wigner distributions as a function of the temperature. 
The relevant parameters are: 7 = 0.25 a;, F — 0.001 lj, ub = 0,0.75,1.5 
respectively represented by dots, circles and asterisks. Full lines are the 
analytical results. 
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5.1.2. Current Since the effect of the oscUlator degree of freedom is entirely included 
in the renormalized rates, the system can be treated formally as a static quantum dot. 
The time-dependent currents thus read: 

(lei 

hit) ^TLPooit) 
In the stationary limit they coincide: 

^stat ^ p^^stat ^ p^^stat ^ bT L ^^g^ 

We show in Figure El the current calculated numerically and the asymptotic value of 
the tunneling regime given by Eq. (|19|l . 

5.1.3. Current-noise We start the calculation with the MacDonald formula for the 
zero frequency current noise [5l ll4[[T5] : 

, oo oo 

S{0) = \im - [ J2 n'Pnit) - ( ^ nP„(t))'] (20) 

ri=0 Ji=0 

where P„(t) is the probability that n electrons have been collected at time t in the 

(n) 

right lead. This probability is connected to the n-resolved probabilities p of the two 
effective states of the tunneling model by the relation: 

Pn{t)=pio\t)+p[f{t) (21) 

(n) 

The n-resolved probabilities pl^ satisfy the equation of motion: 




(22) 



that can be derived by tracing the equation of motion for the total density matrix 
p over bath states with a fixed number (n) of electrons collected in the right lead 
and finally integrating over the mechanical degrees of freedom. The evaluation of 
the different terms of the current-noise (|20|l can be carried out by introducing the 
generating functions Fii{t;z) = ^nPi^\t)z'^ [3 The Fano factor is calculated 
in terms of the stationary probabilities pf,*^* and the pseudoinverse of the Liouvillean 

For a detailed evaluation of the formula H23|) we refer the reader to the section IV B 
of JHI- The resulting Fano factor 

p2 I p2 

F = !^+!« (24) 

assumes the familiar form for a tunneling junction, albeit with renormalized rates. 
In Fig. the value of the Fano factor given by the above formula is depicted as the 
high-damping asymptote of the full calculation. 



Simple models for the quantum shuttle 11 



Current Fano Factor 




0.25 



Figure 6. Left panel - Current as a function of the damping for the 
SDQS. The asymptotic tunneling limit is indicated. The parameters are: 
Tl = Fr = 0.01 oj, 7 = 0.25 a;, d = 0.5 xo, I = 2xo, T = 0. Right 
panel - Current-noise as a function of the damping for the SDQS. The 
asymptotic tunneling limit is indicated. The parameters are the same as 
the ones reported for the current. 



5.2. Shuttling: a classical transport regime 

The simplified model for the shuttling dynamics is based on the observation ~ extracted 
from the full description - that the system exhibits in this operating regime extremely 
low Fano factors {F « 10~^): we assume that there is no noise at all in the system. Its 
state is represented by a point that moves on a trajectory in the device phase-space 
spanned by position, momentum and charge of the oscillating dot. The charge on 
the oscillating dot is a stochastic variable governed by tunnelling processes, however 
in the shuttling regime the tunnelling events are effectively deterministic since they 
are highly probable only at specific times (or positions) defined by the mechanical 
dynamics. 



5.2.1. Equation of motion for the relevant variables We implement the zero noise 
assumption in the set of coupled Klein-Kramers equations (|5J in two steps: we first 
set T = and then simplify the equations further by neglecting all the terms of the 
h expansion since we assume the classical action of the oscillator to be much larger 
than the Planck constant. We obtain: 



dr 



X 



d 



dr 



dP 

Tl _ 
— e 

LO 

X - 



P 



d 

dX 



7 d 
' Z^dP^ 



w; 



cl 
00 



d_ 

dP 



-P 



d_ 

dX 



• dP 



P 



(25) 



ut, X 



(26) 



where we have introduced the dimensionless variables: 

^, P- ^ 
I ' mujl 

The superscript "cl" indicates that we are dealing with the classical limit of the Wigner 
function because of the complete elimination of the quantum "diffusive" terms from 
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the Klein-Kramers equations. In this spirit, it is natural to try an Ansatz for the 
Wigner functions, in which the position and momentum dependencies are separable: 

Wf,iX,P,T)=p,^iT)SiX-X^\T))5iP-P^\T)) 

where the trace over the system phase-space sets the constraint poo + Pii = 1- The 
variables X'^^ and P'^' represent the position and momentum of the (center of mass) 
of the oscillating dot; pii(oo) is the probability for the quantum dot to be charged 
(empty). 

By inserting the Ansatz (|27|) into equation H25|) and matching the coefficients 
of the terms proportional to S x S we obtain the equations of motion for the charge 
probabilities pu: 

Poo = e Poo -i e pn 

^ (28) 
Pu = — e Poo e pii 

LU UJ 

Matching the coefficients proportional to the distributions 6 x S' (here 5' is the 
derivative of the delta-function) yields the equations for the mechanical degrees of 
freedom: 



7 pel) (29) 

LU 

A LU 



(30) 

combined with the normalization condition poo +Pii = 1- Under these conditions the 
system H29|l is equivalent to 

^cl ^ pel 

d -Y (31) 

P^' = -X^' + ^pn - ^P'' 
A uu 

The condition (|30|l also follows by substituting the Ansatz (|27|) into the equations (|25|l 
and by using the equations of motion H31|l and (|28|) . This shows that (|30|l also sets 
the limits of the validity of the Ansatz H27() . However, the only differentiable solution 
for H3()|l is Poo = or pn = for all times, which is not compatible with the equation 
of motion (j2HJ- Thus, the Ansatz 1)27(1 does not yield exact solutions to the original 
equations 1)25(1 . 

While an exact solution has not been found, we can still proceed with the following 
physical argument. Suppose now that the switching time between the two allowed 
states, pii — 1; poo = or poo = 1; Pii — 0, is much shorter than the shortest 
mechanical time (the oscillator period T — 2tt/uj). A solution of the system of 
equations ((31(1 and ((28(1 with this time scale separation would satisfy the condition 
((3U(I "almost everywhere" and, when inserted into ((27(1 would represent a solution for 
(|23. 
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We rewrite the set of equations and as: 
X = P 

P = -X + d*Q --i*P (32) 
Q =rie-2^(l-Q)-r|je2Xg 

where we have dropped the "cl" superscript, renamed pn = Q, used the trace 
condition poo = 1 ^ Pii, and defined the rescaled parameters: d* — d/l, 7* = 
'j/uj, ji — Tl^i^/lo. In the following section we analyze the dynamics implied by 
Eq. (EH). ' 

5.2.2. Stable limit cycles Here we give the results of a numerical solution of Eq. (|32|l 
for different values of the parameters and different initial conditions. For the parameter 
values that correspond to the fully developed shuttling regime, the system has a limit 
cycle solution with the desirable time scale separation we discussed in the previous 
section. Figure [7| shows the typical appearance of the hmit cycle. 

1.5 
1 

S0.5 



810 820 830 -5 5 

cot ^ 

1.5 ■> •= 



O 0.5 



-0.5 

■5 5 -5 5 

X P 

Figure 7. Different representations of the limit cycle solution of the system 
of differential equations i3'^ that describes the shuttling regime. For a 
detailed description see the text. X is the coordinate in units of the 
tunneling length I, while P is the momentum in units of mwl. 

In Fig. [3a) we show the charge Q{t) as a function of time. The charge 
value is jumping periodically from to 1 and back with a period equal to the 
mechanical period. The transition itself is almost instantaneous. In panels (b), (c) 
and (d) three different projections of the 3D-phase-space trajectory are reported and 
the time evolution along them is intended clockwise. The X, P projection shows 
the characteristic circular trajectory of harmonic oscillations. In the X, Q {P, Q) 
projection the position(momentum)-charge correlation is visible. 

The full description of the SDQS in the shuttling regime has a phase space 
visualization in terms of a ring shaped stationary total Wigner distribution function, 
see Fig. |5J We can interpret this fuzzy ring as the probability distribution obtained 
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Figure 8. Correspondence between the Wigner function representation and 
the simplified trajectory limit for the shuttling regime. The white ring is 
the (X, P) projection of the limit cycle. The Q — 1 and Q = portions of 
the trajectory are visible in the charged and empty dot graphs respectively. 
The parameters are 7 = 0.02 oj, d — 0.5 a;o, F — 0.05 lj, I — xq in the upper 
row, 7 = 0.02 cli, d — 0.5 a;o, F = 0.05 lj, I = 2xo in the central row and 
7 = 0.02 Lu, d — 0.5 xo, r — 0.01 uj , I — 2 xo in the lower row. 



from many different noisy realizations of (quasi) limit cycles. The stationary solution 
for the Wigner distribution is the result of a diffusive dynamics on an effective 
"Mexican hat" potential that involves both amplitude and phase of the oscillations. 
In the noise-free semiclassical approximation we turn off the diffusive processes and 
the point-like state describes in the shuttling regime a single trajectory with a definite 
constant amplitude and periodic phase. We expect this trajectory to be the average 
of the noisy trajectories represented by the Wigner distribution. In the third column 
of Fig. IHl the total Wigner function corresponding to different parameter realization 
of the shuttling regime is presented. The white circle is the semiclassical trajectory. 
In the first two columns the asymmetric sharing of the ring between the charged and 
empty states is also compared with the corresponding Q = I and Q = portions of 
the semiclassical trajectory. 

In the semiclassical description we also have direct access to the current as a 
function of the time. For example the right lead current reads: 

lR{r) - Q(r)rfle2^(-) (33) 

and is also a periodic function with peaks in correspondence to the unloading processes. 
The integral of Ir(t) over one mechanical period is 1 and represents the number of 
electron shuttled per cycle by the oscillating dot, in complete agreement with the full 
description. 
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5.3. Coexistence: a dichotomous process 

The longest time-scale in the coexistence regime corresponds to infrequent switching 
between the shutthng and the tunneUng regime, see Fig. El The amplitude of the dot 
oscillations is the relevant variable that is recording this slow dynamics. We analyze 
this particular operating regime of the SDQS in four steps, (i) We first explore the 
consequences of the slow switching in terms of current and current noise, (ii) Next, we 
derive the effective bistable potential which controls the dynamics of the oscillation 
amplitude, (iii) We then apply Kramers' theory for escape rates to this effective 
potential and calculate the switching rates between the two amplitude metastable 
states corresponding to the local minima of the potential, (iv) We conclude the section 
by comparing the (semi) analytical results of the simplified model with the numerical 
calculations for the full model. 



time 



Figure 9. Schematic representation of the time evolution of the current in 
the dichotomous process between current modes in the SDQS coexistence 
regime. The relevant currents are the shuttling (Ish) and the tunneling 
(Itun) currents, respectively. The switching rates Fin and Tout correspond 
to switching in and out of the tunneling mode. 



5.3.1. Two current modes Let us consider a bistable system with two different modes 
that we call for convenience Shuttling (sh) and Tunneling (tun) and two different 
currents /gh and /tun, respectively, associated with these modes. The system can 
switch between the shuttling and the tunneling mode randomly, but with definite 
rates: namely Fin for the process "shuttling — > tunneling" and Fout for the opposite, 
"tunneling — > shuttling" . We collect this information in the master equation: 

For such a system the average current and the Fano factor read |14l [T7| : 

rstat -^shTout + /tunTin 



F = 



Tin Tout 
5(0) _ ^ (/sh - /tun)' TinF, 



(35) 



/^'^* /shTout + /t unTin (Tin + Fout)' 



The framework of the simplified model for the coexistence regime is given by these 
formulas. The task is now to identify the two modes in the dynamics of the shuttle 
device and, above all, calculate the switching rates. This can be done by using the 
Kramers' escape rates for a bistable effective potential. 
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Figure 10. Schematic representation of the effective potentials for the three 
operating regimes. 



5.3.2. Effective potential The tunneling to shuttling crossover visualized by the total 
Wigner function distribution (Fig. [SJl can be understood in terms of an effective 
stationary potential in the phase space generated by the non-linear dynamics of the 
shuttle device. We show in Fig. ^| the three qualitatively different shapes of the 
potential surmised from the observation of the stationary Wigner functions associated 
with the three operating regimes. Recently Fedorets et al. |Sj initiated the study of 
the tunneling-shuttling transition in terms of an effective radial potential. Taking 
inspiration from their work we extend the analysis to the slowest dynamics in the 
device and use quantitatively the idea of the effective potential for the description of 
the coexistence regime. 

In the process of elimination of the fast variables we start from the Klein-Kramers 
equations for the SDQS that we rewrite symmetrized by shifting the coordinates origin 
to d/2: 



dW< 



00 



dt 



dWii 
dt 



d p d 



d 



1 



92 

9p2 



w. 



00 



+ TRe''^/'W,^-TLe-^^/'Y.T^^ 
d 



I J 9p2n 



dp^ 



(36) 



Wii 



V (-1)" 
^ (2n)! 

n=0 ^ ' 



where the renormalization of the tunneling rates due to the coordinate shift has been 
absorbed in a redefinition of the F's. The idea is to get rid of the variables that due 
to their fast dynamics are not relevant for the description of the coexistence regime. 
In equations 1)36(1 we describe the electrical state of the dot as empty or charged. We 
switch to a new set of variables with the definition: 



W± = Woo ± Wii 



(37) 



In absence of the harmonic oscillator the state |-|-) would be fixed by the trace sum rule 
and the state |— ) would relax to zero on a time scale fixed by the tunneling rates. We 
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assume that also in the presence of the mechanical degree of freedom the relaxation 
dynamics of the |— ) state is much faster than the one of the |+) state. 

In the dimensionless phase space given by the coordinates X and P of H2t)|) we 
switch to the polar coordinates defined by the relations [S]: 

X ^ Asincj) P = Acos(l) (38) 

Since we are interested only in the dynamics of the amplitude in the phase space (the 
slowest in the coexistence regime) we introduce the projector that averages over 
the phase : 

r^']^^ f d^' (39) 

We also need the orthogonal complement = 1 — 7^0. Using these two operators we 
decompose the Wigner distribution function into: 

W+ = r4,W+ + Q^W+ ^W+ + W+ (40) 

Finally we make a perturbation expansion of 136|l in the small parameters: 

These three inequalities correspond to the three physical assumptions: 

(i) The external electrostatic force is a small perturbation of the harmonic oscillator 
restoring force in terms of the sensitivity to displacement of the tunneling rates. 
This justifies an oscillator-independent treatment of the tunneling regime. 

(ii) The tunneling length is large compared to the zero point fluctuations. Since 
the oscillator dynamics for the shuttling regime (and then partially also for the 
coexistence regime) happens on the scale of the tunneling length, this condition 
ensures a quasi-classical behaviour of the harmonic oscillator. 

(iii) The coupling of the oscillator to the thermal bath is weak and the oscillator 
dynamics is under-damped. 

Using these approximations the Klein-Kramers equations (|36|l reduce (for details 
see, e.g., [Sl^J) to the form: 

drW+{A,T) = jdAA[V'{A) + D{A)dA]W+{A,T) (42) 

where V'{A) = ^V{A) and D{A) are given functions of A. Before calculating 
explicitly the functions V and D we explore the consequences of the formulation 
of the Klein-Kramers equations H3ti|l in the form H42|) . The stationary solution of the 
equation reads |S]: 

»-..(^)4exp(-^^^'^] (43) 

where Z is the normalization that ensures the integral of the phase-space distribution 
to be unity: dA'2T:A'W^''\A') = 1. Equation ^ is identical to the Fokker- 
Planck equation for a particle in the two-dimensional rotationally invariant potential 
V (see Fig. I1(J|I with stochastic forces described by the (position dependent) diffusion 
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coefficient D. All contributions to the effective potential V and diffusion coefficient D 
can be grouped according to the power of the small parameters that they contain. 

2 



«2(A) + ^^a3(A) 
uj 21 



(44) 



where the a functions read: 

ao —T-'4> cos (/)Gor_ 
1 



ai = - -V^cos(t)T^dp{GoT^) 
Q!2 =:P0Cos(?!)Gor_goQ09p(Gor_) 



(45) 



"3 COS 

and the /3's can be written as: 
Pi =^7^0 cos^ ( 



A . 

— sm ( 
2 



'(Gor_) 



r+-r_Gor- 



P2 ^V^icoscf) Gocos0 + Gor_goQ0Cos0Gor_ 
/33 ^AV4,cos(t)\Go cos^ (/)Gor_ + ^Gor^ sin20 ~ isin20Gor 



(46) 



where 



r+ - Tre 



90 ^{d^r' (47) 

Go = (90 + r+)-i 

The a and (3 functions are calculated by isolating in the Liouvillean for the distribution 
W+ the driving and diffusive components with generic forms 



jdAA{aM)} 



and 



A 



dAA{p,{A)}dA 



(48) 



(49) 



respectively. As an example we give the derivation of the functions 013 and fi^. We 
start by rewriting the equation of motion H42|) for the distribution Wl|_ in the form: 

drW+ = £[W+] « {C^ + C^^)[W+] (50) 

where we have distinguished the Liouvilleans of first and second order in the small 
parameter expansion (|41l) . The contribution ^ f of the second order Liouvillean 
reads: 



C^d = V^dp[GodpPGoT^ + PgoQ4,dpGor_ + GoT^goQ^dpP] 



(51) 
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and represents the starting point for the calculation of the functions and (i^. We 
express then the differential operators dp in polar coordinates and take into account 
that the Liouvillean is applied to a function independent of the variable (f): 



£^rf ^^dAAV^coa(j)\GlT^ + GoAcos(j)dp{GoT-) + GoAcos^ cjyGoT-dA 

+ ^cos0goQ0C^p(Gor-) + Acos05oQ0 cos0Gor„9yi 
+ Gor-goQ0^cos2 <\>d 



(52) 



Finally, we separate in (|52|l the driving term from the diffusive contributions and thus 
identify the functions and %: 



1 

A' 



OaA < Vd, cos ( 



Glr^+AGodp{GoT^) ^ ^ sm<l>dp{Gor^] 



Go cos^ 0Gor_ + ^GqT^ sin 20 - ^ sin 2(f>GoT_ 



(53) 



dA 



Some of these results appear also in the work by Fedorets et al. 5^. Since we have 
projected out the phase we are effectively working in a one-dimensional phase 
space given by the amplitude A. Note however that Eq. H42|l though is not as it 
stands a Kramers equation for a single variable. This is related to the fact that also 
the distribution W+ is not the amplitude distribution function, but, so to speak, 
a cut at fixed phase of a two dimensional rotationally invariant distribution. The 
difference is a geometrical factor A. We define the amplitude probability distribution 
W(A, r) — AW+{A,t) and insert this definition in equation H42() . We obtain: 

drW{A, t) = dAA[V'{A) + D{A)dA] jW{A, r) 
^ dA[V{A) + D{A)dA]WiA,T) 
where we have defined the geometrically corrected potential 

'^DiA') 



V(A) = V{A) 



Ao 



A' 



-dA' 



(54) 



(55) 



which for an amplitude independent diffusion coefficient gives a corrected potential 
diverging logarithmically in the origin. The lower limit of integration is arbitrary and 
reflects the arbitrary constant in the definition of the potential. The equation (|54|l 
is the one-dimensional Kramers equation that constitutes the starting point for the 
calculation of the switching rates that characterize the coexistence regime. 

The effective potential V that we obtained has, for parameters that correspond 
to the coexistence regime, a typical double- well shape (see e.g. the left panel of figure 
lll|) . We assume for a while the diffusion constant to be independent of the amplitude 
A. In this approximation the stationary solution of the equation (|54|l reads: 



1 

z 



cxp 



V{A) \ 
D / 



(56) 



where Z is the normalization: Z = exp ^— ^ dA. The probability distribution 
is concentrated around the minima of the potential and has a minimum at the potential 

J In this step of the derivation we have also used the projector 7-"^ to define a scalar product 
'P.pfWgW = (/iff) and the adjoint relation: {f,dg) = {O^ f,g). 
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Figure 11. Left panel - Bistable effective potential for the SDQS coexistence 
regime. The important amplitudes for the calculation of the rates are 
indicated. In red (blue) are indicated the reflecting (full) and absorbing 
(dashed) borders for the calculation of the Fout ^inj escape rate. Right panel 
- Example of the stationary distribution W^^^ (blue) and the amplitude 
distribution W"*''' (green) for the SDQS in the coexistence regime. The 
tunneling and shuttling states are in both cases well separated. 



barrier (see e.g the right panel of Fig. If this potential barrier is high enough 

{i.e. Vmax ~ Vniin 3> D) wB clearly identify two distinct states with definite average 
amplitude: the lower amplitude state corresponding to the tunneling regime and the 
higher to the shuttling. 

The coexistence regime of a SDQS is mapped into a classical model for a particle 
moving in a bistable potential V with random forces described by the diffusion constant 
D. The correspondent escape rates from the tunneling to the shuttling mode (Fout) 
and back (Fin) can be calculated using the standard theory of Mean First Passage 
Time (MFPT) for a random variable (18j : 




-1 

dBe"^ I dAe^"^ 



dBe o t dAe o 
J B 



(57) 



where integration limits of equation (|57|) are graphically represented in the left panel 
of Fig. 111! We can now insert the explicit expression for the switching rates Fin and 
Fout in Eq- H35|) and obtain in this way the current and Fano factor for the coexistence 
regime. They represent, together with the stationary distribution H5()|l the main result 
of this section and allow us for a quantitative comparison between the simplified model 
and the full description of the coexistence regime. 



5.3.3. Comparison The phase space distribution is the most sensitive object to 
compare the model and the full description. One of the basic procedures adopted 
in the derivation of the Kramers equation (|54|l is the expansion to second order in 
the small parameters H41II . In order to test the reliability of the model we simplify 
as much as possible the description reducing the model to a classical description: 
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namely taking the zero limit for the parameter (^)- We reahze physicahy this 
condition assuming a large temperature and a tunneling length / of the order of the 

thermal length Ith — \J~^^- Also the full description is slightly changed, but not 
qualitatively: the three regimes are still clearly present with their characteristics. The 
numerical calculation of the stationary density matrix is though based on a totally 
different approach. In the quantum regime we used the Arnoldi iteration scheme 
for the numerically demanding calculation of the null vector of the big (10^ x 10"*) 
matrix representing the Liouvillean 01 Ej. Problems concerning the convergence of 
the Arnoldi iteration due to the delicate issue of preconditioning forced us to abandon 
this method in the classical case. We adopted instead the continued fraction method 
[TS). In Figs. lT^IT^ andiniwe present the results for the stationary Wigner function, 




Figure 12. Stationary amplitude probability distribution W for the SDQS in 
the coexistence regime. We compare the results obtained from the simplified 
model (full line) and from the full description (asterisks). These results are 
obtained in the classical high temperature regime ksT ^ ftw. The amplitude 

is measured in units of /th = \J • The mechanical damping 7 in units 
of the mechanical frequency u). The other parameter values are d = 0.05 /th 
and r = 0.015 LJ, A = 2Ath. 



the current and the Fano factor, respectively, in the semiclassical approximation and 
in the full description. 

Deep in the quantum regime the coexistence regime {e.g. Fig. [21 where the 
amplitude of the shuttling oscillations is « 7xo) is not captured quantitatively with 
the simplified model. Given that the concept of elimination of the fast dynamics 
is still valid, we believe that the discrepancy indicates that the expansion in the 
small parameters has not been carried out to sufficiently high order . The effective 
potential calculated from a second order expansion still gives the position of the 
ring structure with reasonable accuracy but the overall stationary Wigner function 
is not fully reproduced due to an inaccurate diffusion function D{A). One should 
thus consider higher order terms in the parameter (xo/l)^. A higher order expansion, 
however, represents a fundamental problem since it would produce terms with higher 
order derivatives with respect to the amplitude A in the Fokker-Planck equation and, 
consequently, a straight forward application of the escape time theory is no longer 
possible. 

It has nevertheless been demonstrated |9j with the help of the higher cumulants 
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Figure 13. Current in the coexistence regime of SDQS. Comparison between 
semianalytical and full numerical description. For the parameter values see 



• numerical 
semianalytical 





Figure 14. Fano factor in the coexistence regime of SDQS. Comparison 
between semianalytical and full numerical description. For the parameter 
values see Fig. 



of the current that the description of the coexistence regime as a dichotomous process 
is vahd also deep in the quantum regime {I = 1.5 xq), the only necessary condition 
being a separation of the ring and dot structures in the stationary Wigner function 
distribution. 

We observe that the second order expansion for the effective potential H41|) is 
essentially converged, and able to give the correct position of the shuttling ring also 
in the quantum regime. From Eq. (|43|l it is clear that a strongly amplitude dependent 
diffusion constant D{A) would destroy this agreement. We conjecture that a higher 
order expansion may lead to an effective renormalization of the diffusion constant. 
To test this idea we used the diffusion constant as a fitting parameter, and found 
that the current and the Fano factor are very accurately reproduced by using a fitted 
diffusion constant, with a value approximately twice larger than the one calculated at 
zero temperature and in second order in the small parameters. Clearly, further work 
is needed to find out whether the agreement is due to fortuitous coincidence, or if a 
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physical justification can be given to this conjecture. Also, an investigation of whether 

the results obtained in the classical limit are extendable to the quantum limit by a 
controlled renormalization of the diffusion constant is called for. 

6. Conclusions 

The specific separation of time scales in the diS'erent regimes allowed us to identify 
the relevant variables and describe each regime by a specific simplified model. In the 
tunneling (high damping) regime the mechanical degree of freedom is almost frozen and 
all the features revealed by the Wigner distribution, the current and the current noise 
can be reproduced with a resonant tunneling model with tunneling rates renormalized 
due to the movable quantum dot. Most of the features of the shuttling regime (self- 
sustained oscillations, charge-position correlation) are captured by a simple model 
derived as the zero-noise limit of the full description. Finally for the coexistence regime 
we proposed a dynamical picture in terms of slow dichotomous switching between 
the tunneling and shuttling modes. This interpretation was mostly suggested by the 
presence in the stationary Wigner function distributions of both the characteristic 
features of the tunneling and shuttling dynamics and by a corresponding gigantic 
peak in the Fano factor. We based the derivation of the simplified model on the fast 
variables elimination from the Klein-Kramers equations for the Wigner function and 
a subsequent derivation of an efi'ective bistable potential for the amplitude of the dot 
oscillation (the relevant slow variable in this regime). 

The comparison of the results obtained using the simplified models with the full 
description in terms of Wigner distributions, current and current-noise proves that 
the models, at least in the limits set by the chosen investigation tools, capture the 
relevant features of the shuttle dynamics. 

The work of T. N. is a part of the research plan MSM 0021620834 that is financed 
by the Ministry of Education of the Czech Republic, and A. D. acknowledges financial 
support from the Deutsche Forschungsgemeinschaft within the framework of the 
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